1 はじめに

RをGISとして使う. ポリゴンオブジェクト上に位置するラインオブジェクトの延長を算出する. それぞれ以下のようにする.

  • ポリゴン:土砂災害警戒区域
  • ライン:道路

土砂災害警戒区域上に位置する道路を特定し,その延長を計算することで防災政策を立案するうえでの一つの判断材料を与える.

2 前準備

使用するライブラリの読み込み,国土数値情報からのダウンロード,openstreetmapからのダウンロードを行う.

2.1 ライブラリの読み込み

library(sf)
library(tidyverse)
library(tmap)
library(mapview)

2.2 国土数値情報から

石川県のやつについて,以下を読み込む.

  • 行政区域:N03~
  • 緊急輸送道路:N10-15~
  • 土砂災害警戒区域:A33-17~

コードは以下のとおりである.

# 日本語を扱うのでCP932でのエンコードを忘れない
isikawa <- st_read("N03-170101_17_GML/N03-17_17_170101.shp", options = "ENCODING=CP932", stringsAsFactors=F)
emergency <- st_read("N10-15_17_GML/N10-15_17.shp", options = "ENCODING=CP932", stringsAsFactors=F)
sediment <- st_read("A33-17_17_GML/A33-17_17_GML/A33-17_17Polygon.shp", options = "ENCODING=CP932", stringsAsFactors=F)

N03_004が市町村名に相当する. ここを七尾市に指定すればよさそうである.

nanao <- isikawa %>%
  filter(N03_004 == "七尾市")
## Warning: package 'bindrcpp' was built under R version 3.4.4

行政区域を七尾市のみに絞ったものができた. このポリゴン上で,道路とか土砂災害警戒区域を絞り込む.

emergency_nanao <- st_intersection(nanao, emergency)
sediment_nanao <- st_intersection(nanao, sediment)

そして投影座標系をEPSG2449とする.

nanao <- st_transform(nanao, 2449)
emergency_nanao <- st_transform(emergency_nanao, 2449)
sediment_nanao <- st_transform(sediment_nanao, 2449)

mapview()を使ってインタラクティブに確認する. 車線名(国道n号など)に応じて道路の色を変更する.

mapview(nanao) +
  mapview(emergency_nanao,
          zcol = "N10_004")
#colnames(emergency_nanao)

2.3 openstreetmapから

七尾市周辺のbboxについて,以下を読み込む.

  • morotway
  • trunk
  • primary
  • secondary
  • tertialy
  • unclassified
  • residential

openstreetmapのapiをつかう. pythonのラッパーがあるので,それを使用する. Rからpythonを呼び出すにはこのパッケージをつかう.

bboxの範囲は(36.964677, 136.772596,37.198499, 137.070257)である. pythonを使う下準備は次の通りとなる.

reticulate::use_python(Sys.which(names = "python"))
opass <- reticulate::import(module = "overpass")
api = opass$API() # Rからは$で潜っていく.

たくさん作るので,関数にしておく

geojson_as_sf <- function(way_attr){
  # way_attr は文字列
  bound_box <- "36.964677, 136.772596,37.198499, 137.070257"
  query_opass <- str_c('way["highway"="',way_attr,'"](',bound_box,');(._;>;)')
  attr_api <- api$get(query_opass, verbosity = 'geom')
  way_sf_NodeWay <- st_as_sf(st_read(attr_api))
  way_sf <- way_sf_NodeWay %>%
    filter(highway == way_attr)
  way_sf <- st_transform(way_sf, 2449) # 4326から2449にする
}

道路属性ごとに異なるsfオブジェクトを生成する. 道路属性を文字列として格納しgeojson -> sf する.

# やたらと出力が長くなる. message = F だとhtml出力が止まる.
attr_str <- c("motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential")
for (i in attr_str) {
  now_sf <- str_c(i, "_sf")
  assign(now_sf, geojson_as_sf(i))
}

バウンディボックスは矩形であり,七尾市の外側の道路も取得している. これを七尾市上のみに絞り込む.

for (i in attr_str) {
  now_sf <- str_c(i, "_sf")
  nanao_sf <- str_c(i, "_nanao") # 道路のsfについて 七尾市だけのものをこれに格納する
  assign(nanao_sf, st_intersection(nanao, get(now_sf))) # nanao_sf <- st_intersection(nanao, atr_sf) とするのをまあ
}

最終的に,全道路属性が1つのデータフレームにまとまっているsfオブジェクトを生成する. 現時点ではまだ,道路属性ごとに異なるsfオブジェクトである. これを以下の手順でひとまとめにする.

  1. 道路属性個数分だけあるsfオブジェクトをデータフレームとsfcオブジェクトに分解.
  2. 1.で,道路属性個数分のデータフレームができたので,それをひとまとめにする.
  3. 1.で,道路属性個数分のsfcオブジェクトができたので,それをひとまとめにする.
  4. sfをつくる.2.でできた1つのデータフレームと,3.でできた1つのsfcをくっつける

分解作業,すぐおわる.

for(i in attr_str){
  now_nanao <- str_c(i, "_nanao")
  
  # データフレームだけを now_df に格納
  now_df <- str_c(i, "_df")
  assign(now_df, st_set_geometry(get(now_nanao), NULL))
  
  # sfc だけを now_sfc に格納
  now_sfc <- str_c(i, "_sfc")
  assign(now_sfc, st_geometry(get(now_nanao)))
}

bind_rowにより,道路属性ごとのデータフレームを下に付け足す.

# 格納用のデータフレームを作って,それにバインドしていく
Road_attr_df <- data.frame()
for(i in attr_str){
  now_df <- get(str_c(i, "_df")) # ここでget()すれば1行でデータフレームになる
  Road_attr_df <- bind_rows(Road_attr_df, now_df)
}

sfcも全道路属性をひとまとめにする. リストに格納する.

# 空のリストに放り込む
Road_sfc <- list()

for (i in attr_str) {
  now_sfc <- get(str_c(i, "_sfc")) # ここでget()すれば1行で sfc
  Road_sfc <- c(Road_sfc, i = now_sfc)
  
}

sfオブジェクトをつくる.

Road_sf <- st_sf(Road_attr_df, geometry = Road_sfc)
Road_sf <- st_set_crs(Road_sf, 2449) # 本当は上のやつでいっしょにしたい

道路属性ごとに色分けして図化する.

mapview(nanao) +
  mapview(Road_sf,
          zcol = "highway")

3 土砂災害警戒区域上に限定する

highway属性をすべて有する道路オブジェクトRoad_sfが作られた. 土砂災害警戒区域の領域に含まれる道路延長の算出のため,Road_sfをこの領域内だけに限定したオブジェクトをさらに作成する. 土砂災害警戒区域のsfオブジェクトはsediment_nanaoである.

Road_sedi <- st_intersection(sediment_nanao, Road_sf)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

bboxをいきなりでなく,sfcでよい??

p1 <- st_point(c(136.951, 37.001))
p2 <- st_point(c(136.980, 37.025))
interest_sfc <- st_sfc(p1, p2, crs = 4326)
interest_sfc <- st_transform(interest_sfc, 2449)
interest_sfc
## Geometry set for 2 features 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -19194.69 ymin: 111089.9 xmax: -16608.42 ymax: 113747.7
## epsg (SRID):    2449
## proj4string:    +proj=tmerc +lat_0=36 +lon_0=137.1666666666667 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## POINT (-19194.69 111089.9)
## POINT (-16608.42 113747.7)

特定の範囲のみを拡大して表示する.

tm_shape(nanao, bbox = interest_sfc) +
  tm_polygons() +
tm_shape(sediment_nanao) +
  tm_polygons(col = "red") +
tm_shape(Road_sf) +
  tm_lines() +
tm_shape(Road_sedi) +
  tm_lines(lwd = 5)

4 距離を算出する

ラインオブジェクトの長さの算出はst_length()で実施する.

# 全道路延長の算出
Road_sf <- Road_sf %>%
  mutate(road_length = st_length(Road_sf))
# highway属性ごとに集約して算出
Road_sf_sum <- st_set_geometry(Road_sf, NULL) %>%
  group_by(highway) %>%
  summarise(ALL = sum(road_length))

# 土砂災害警戒区域上の道路延長の算出
Road_sedi <- Road_sedi %>%
  mutate(road_length = st_length(Road_sedi))
# highway属性ごとに集約して算出
Road_sedi_sum <- st_set_geometry(Road_sedi, NULL) %>%
  group_by(highway) %>%
  summarise(sedi = sum(road_length))

土砂災害区域に含まれる道路延長は全道路延長の何パーセントに相当するか.

Road_sedi_sum$sedi / Road_sf_sum$ALL * 100
## Units: [1]
## [1]  4.510644  5.863741  8.557915 10.498615  7.238726  8.416303  9.448384

表形式で整理する.

full_join(Road_sedi_sum, Road_sf_sum) %>%
  mutate(ratio = round(sedi/ALL*100,2)) %>%
  knitr::kable()
## Joining, by = "highway"
highway sedi ALL ratio
motorway 1790.554 [m] 39696.19 [m] 4.51 [1]
primary 3166.712 [m] 54004.97 [m] 5.86 [1]
residential 43276.691 [m] 505692.01 [m] 8.56 [1]
secondary 15148.003 [m] 144285.73 [m] 10.50 [1]
tertiary 10674.540 [m] 147464.35 [m] 7.24 [1]
trunk 4234.333 [m] 50311.08 [m] 8.42 [1]
unclassified 33719.617 [m] 356882.39 [m] 9.45 [1]